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ABSTRACT 

Context. The non-linear Hall term present in the induction equation in the electron-magneto-hydrodynamics limit is responsible for 
the Hall drift of the magnetic field and, in some cases, for the so-called Hall instability. 

Aims. We investigate whether or not the growth rates and eigenfunctions found in the linear analysis are consistent with the results of 
non-linear numerical simulations. 

Methods. Following the linear analysis of Rheinhardt & Geppert, we study the same cases for which the Hall instability was predicted 
by solving the non-linear Hall induction equation using a two-dimensional conservative and divergence-free finite difference scheme 
that overcomes intrinsic difficulties of pseudo-spectral methods and can describe situations with arbitrarily high magnetic Reynolds 
numbers. 

Results. We show that unstable modes can grow to the level of the background field without being overwhelmed by the Hall cascade, 
and cause a complete rearrangement of the field geometry. We confirm both the growth rates and eigenfunctions found in the linearized 
analysis and hence the instability. In the non-linear regime, after the unstable modes grow to the background level, the naturally 
selected modes become stable and oscillatory. Later on, the evolution tends to select the modes with the longest possible wavelengths, 
but this process occurs on the magnetic dilfusion timescale. 

Conclusions. We confirm the existence of the Hall instability. We argue against using the misleading terminology that associates the 
non-linear Hall term with a turbulent Hall cascade, since small-scale structures are not created everywhere. The field evolves instead 
in a Burgers-like manner, forms local structures with strong gradients which become shocks in the zero resistivity limit, and Hall 
waves are launched and propagated through the entire domain. 
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1. Introduction 

The influence of the Hall term on the magnetic field evolution 
. of a neutron star (NS) has been debated for many years. By 
' analogy with the vorticity equation in ordinary hydrodynamics, 
I iGoldreich & Reise negger (1992) conjectured that the transfer of 

■ magnetic energy from large to small scales proceeds in a simi- 
lar way to ordinary turbulence. However, the analogy of the Hall 
induction equation with the vorticity equation is not complete, 
and the conjecture remained to be confirmed by multidimen- 
sional numerical simulations. V ainshtein et al.i (2000) proposed 
a mechanism for the fast dissipation of magnetic field based on 

■ the Hall drift in stratified media. They correctly pointed out that 
' Hall currents are able to create current sheets (which are sites for 

efficient dissipation) and that the evolution of the toroidal field 
resembles the Burgers equation. The same Burgers-like equation 
is applicable even to no n-stratified media, but in a spherical shell 
(iPons & Geppertll2007l) . in which the Hall term in the induction 
equation tends to create current sheets instead of ordinary tur- 
bulence. Rheinhardt & Geppert (2002) (hencefort RG) showed 
by a linear analysis that, in a one-component (electron) plasma, 
a large-scale background magnetic field may become unstable 
to smaller scale perturbations. This Hall-drift induced instabil- 
ity (HI) occurs when the magnetization parameter is high and 
the background field has enough curvature. Since these condi- 
tions may be realized in the crust of a NS, the problem of the 
HI became interesting not only from the MHD point of view 



but also for the astrophysics community. Although the Hall-drift 
is a non-dissipative process, the growth of small-scale magnetic 
field components modifies the overall magnetic field structure 
and opens the possibility of more rapid field decay than pure 
ohmic dissipation would predict. This effect may be relevant, at 
least, in relatively young neutron stars. Determining its effect on 
the long-term evolution, and which configurations ar e stable on 
long time-scales are different issues. A review by Reis eneggeJ 
(i2009.) discusses the main properties of stable magnetic config- 
urations emphasizing the importance of a non-barotropic strati- 
fication in neutron stars, and provides insight into the expected 
long-lived magnetic field structures. Some observational impli- 
cations of rapid magnetic fi e ld decay were also discussed by 
'Geppert & Rheinhardt! (l2002h : lPons etan(l2007h : I Aguilera et all 
(2008b a). Other attempts to understand the long-term evolu- 
tion of neutron star magnetic fields by numerical simulations 
were performed by Hovo s et al.. (2008) . but were restricted to 
1 -dimensional geometries. 

The appeara nce of the HI was disputed by 
Hollerbach & Riidigeil ( 120021) : they did not find any evi- 
dence of the instability in their numerical study of field decay 
under the influence of the Hall drift in a constant density 
spherical shell. The study was later ext ended to include density 
gradients (iHollerbach & Rudigeill2004l) . However, it was clear 
that they could not find the HI because of numerical limita- 
tions of their code, which is unable to handle situations with 
sufficiently high magnetic Reynolds number. The relevance of 
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the HI to astr ophysical scenar i os had akeady been questioned 
by RG. Later ICumming et alJ (l2004l) conjectured that the gap 
in the energy spectrum created by the HI could be rapidly 
filled in by the Hall cascade. They also computed Hall-wave 
eigenfunctions, including the effect of density gradients across 
the crust, and found agreement with the results of RG. An 
important component of their study was their estimate of the 
respons e of the crust to the ma gnetic stresses induced by Hall 
waves. iPons & GeppertI ( l2007h performed 2D simulations by 
solving the non-linear Hall-induction equation in a realistic 
NS crust and found some evidence for a rapid reorganization 
of the field during a relatively short period, but no clear con- 
firmation of the HI could be concluded. The code again had 
numerical limitations that avoided additional exploration of 
extreme limits to confir m or rule out the existence of the HI. 
Wareing & Hollerbach (l2009h rhereafter WH) used their 3D 
pseudo-spectral code and concluded that the HI proposed by RG 
is a consequence of particular initial conditions and irrelevant 
compared to the Hall cascade. They considered a similar, but not 
identical, geometry as RG and used a code limited to employing 
periodic boundary conditions in all directions. 

In this Letter, we solve the non-linear induction equation us- 
ing a new code (Pons et al., 2010, in preparation) based on a 
divergence-free, conservative, finite-difference scheme that re- 
sembles numerical methods widely employed in astrophysical 
MHD problems to handle situations where shocks are present. It 
solves problems caused by the transition from a parabolic equa- 
tion to the hyperbolic regime as the magnetic resistivity van- 
ishes. Since the code is designed by following methods devel- 
oped to solve hyperbolic equations, which are similar to stan- 
dard te chniques used in p erfect fluid hydrodynamics and ideal 
MHD dAnton et al.ll2006h . it is able to deal with zero physical 
resistivity, therefore extending all previous studies to the limit 
of arbitrary high magnetic Reynolds numbers. Our numerical 
scheme leads to a fourth order numerical hyper-resistivity, in 
a similar way to other s tudies on the basis of which we devel- 
oped our own algorithm dO' Sullivan & Downesl2006l:lT6th et al] 
l2008h . We performed a number of tests including the purely dif- 
fusive case and the exactly non-resistive case (e.g., propagation 
of whistler waves on a homogeneous background). In all cases 
where analytical solutions are available, the correct modes, de- 
cay times, and wave propagation speeds were reproduced. 

This allows us to study the same model originally analyzed 
by RG and answer open questions about the occurrence of the 
HI and its interplay with the so-called Hall cascade. 

2. Basic equations and model set up. 

To be able to compare with the results of RG, we restrict our- 
selves to the case of constant density and electric conductivity. 
In dimensionless variables (for details, see RG), the Hall induc- 
tion equation is given by 

— = AB- Vx[(VxB)xB] , (1) 

ot 

which must be solved with the constraint V ■ B - Q. 

Following RG, we consider a slab that is infinitely extended 
in both X- and y-directions. The upper boundary is located at 
z = +1 in contact to vacuum, while the lower boundary, lo- 
cated at z = -1, forms the interface with a perfect conduc- 
tor. The background (initial) magnetic field is assumed to have 
only the jc-component and depends only on the z-coordinate 
^0 - B()f{z)ex, where e.i is the unit vector in the x-direction. 



Bo is a constant fixing the maximum field strength, and f(z) is 
an arbitrary function. The reason for this choice is that we must 
start with a background field that is not affected by the non-linear 
Hall term. Note that this is not an arbitrary choice: in any prob- 
lem where one wants to study the evolution of perturbations on 
a background, the first important requirement is that the back- 
ground has to be an equilibrium solution. Any other choice of 
background field, will cause rapid dynamical evolution. After 
all, nobody would attempt, for instance, to study perturbations 
within a star starting with a fluid configuration in which gravity 
is not balanced by the gradient of pressure. Thus, imposing the 
condition Vx[(V x Bo) x Bq] - 0, which was criticized by WH, 
is a pre-requisite for a linear analysis, not just a special choice. 

A different question is whether or not this background is 
an accurate choice and how this affects the astrophysical rele- 
vance of the conclusions. Naively, since a NS at birth is fluid, 
one would expect that an MHD equilibrium is reached before 
the crust is formed, and thus that the initial configuration may 
be closer to an MHD equilibrium than to a Hall equilibrium. 
But the conditions for reaching these two equilibria are simi- 
lar (not identical) and one may expect that the initial model is 
not far from a Hall equilibrium. Another astrophy sic ally inter- 
esting scenario could be relatively old NS, with an age of about 
10^ years. During the photon cooling era, the temperature can 
become so low that the magnetization parameter becomes high 
even for moderate magnetic fields (10'^ - lO'-' G), but it does 
take time for the temperature to reach such low values. The sys- 
tem slowly approaches slowly the non-resistive case, for a rel- 
atively long period not being affected by the HI. It is reason- 
able to assume that the crustal magnetic field has enough time 
to steadily adjust itself to a Hall equilibrium until the conditions 
for triggering the instability can be reached in a similar manner 
as described above. 

The choice of the function f(z) is of course not unique. To 
be able to compare with RG, we use 

f(z)^(l+z)(l-zh (2) 

for which the modes with the most rapid growth rate correspond 
to horizontal wavenumbers ky = and ky- of order 1 for Bo ~ 
1000. Therefore, we restrict ourselves to the 2D case with ky = 0, 
e.g., we assume independence on the y-coordinate. 

All quantities are assumed to depend periodically on x with 
period L (that is, the length of the numerical domain in the x- 
direction), and periodic boundary conditions are applied in terms 
of the jc-coordinate. For the inner boundary condition, as in for- 
mer studies, we mimic the conditions in a NS's crust by as- 
suming that the matter beneath the lower surface at z = -1 is 
a perfect conductor. In that case, the Meissner-Ochsenfeld ef- 
fect prevents the magnetic field from penetrating into the super- 
conducting matter. This requires that the normal component of 
the magnetic field and the tangential comp onents of the elec- 
tric field have to vanish (iPons & Gepperfll2007l) . i.e., we assume 
that B- - Ex - Ey = at that interface. At the transition to 
a vacuum (z - +\), the analog to the toroidal component By 
has to vanish, but the boundary conditions for the other compo- 
nents are not so simple. In terms of individual modes with fixed 
wavenumber, a spectral representation that matches the vacuum 
solution needs to satisfy a simple relation between the Stokes 
function and its derivative (RG). However, in real space and us- 
ing the components of the magnetic field, we have to employ a 
non-local boundary condition. By using Green's representation 
formula to solve the Laplace equation in a bounded domain, and 
considering that we reduce the problem to 2D, one can derive an 
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explicit form for the outer boundary condition (Rheinhardt 2009, 
private communication) given by 



- ^^^^^ E — 771 



\k=-N ' 



, at z - +1 



(3) 



where the integer N > must be chosen to be high enough to 
ensure that the sum converges. 

With the same initial model and physically motivated bound- 
ary conditions, RG found that the HI appears for Bq > 3, and 
that some small-scale structures must grow quickly. Since there 
is some confusion in the literature about the interpretation of the 
term small scales, we briefly review their results and make some 
remarks. For sufficiently strong and curved background fields, 
there is always a region in the wavenumber space where unstable 
modes appear. The highest growth rates correspond to ky - and 
the range of allowing for unstable modes changes depending 
on both the form of f(z) and the value of Bo but, in dimension- 
less units and with the slab defined in the region z = [-1,-1-1], 
it is in all cases of order 1, or even smaller. This means that the 
most rapidly growing unstable mode has a typical wavelength of 
the order of the slab height. Thus, in some sense, it cannot be 
rigorously assumed to be small-scale mode. However, we note 
that the typical crust size of a NS is about 1 km, compared to 
its radius of about 10 km. Thus, if the HI occurs, one would ex- 
pect structures of a typical size of one km to grow within the 
crust, with a typical area of a; 10^ km^. In this respect, the field 
geometry will differ considerably from a typical dipolar field, 
consisting instead of tens or hundreds of smaller structures, and 
that is what RG meant by small-scale structures. In addition, 
the toroidal field (here B,,) may grow rapidly in a very narrow 
region of 10-100 m beneath the surface. This is far from the 
typical picture of turbulence in fluid dynamics, where multiple 
small-scale substructures grow in a large region at the same time. 
In the HI, there are a few small-scale structures, which can be 
very localized, in particular close to the boundaries, or wherever 
the toroidal field changes sign, because the Burgers-like term in 
the evolution equation will produce a sharp discontinuity. This 
qualitative difference is easy to visualize when one follows the 
evolution of the field components in the real space. It becomes, 
however, more difficult to distinguish when working in wave- 
number space with pseudo-spectral methods if we only consider 
the power-law in the spectrum. We think that most of the dis- 
agreement and confusion in the past few years has been caused 
by an inaccurate use of language, terms such as Hall cascade or 
small scale being interpreted differently by different readers. 

3. Results 

WH claimed that any instabilities are always overwhelmed by a 
turbulent Hall cascade. To either confirm or reject these conclu- 
sions, we performed simulations following the non-linear evolu- 
tion of the initial model described above. We assumed that Bq - 
1000, which, according to the linear analysis, should be affected 
by the HI. The wavenumber of unstable modes is in the range 
k^ - [0.3, 3] and growth rates in the range 2 x lO""* - 10"^ ohmic 
decay times. We perturbed the background field with small per- 
turbations (0.1-1% in amplitude) of different wavenumbers in 
that range, and followed the time evolution for many growth 
times, until a quasi-stationary diffusion regime was reached. In 
most cases, we found that the growth of the structures closely 
resembled the eigenfunctions obtained in RG with the correct 
timescales. Figure [1] shows the results for an initial perturbation 
of kx - 0.5 and maximum amplitude 1. This mode should be 



unstable and its growth time is about 2 x 10"^. As predicted, 
the initial perturbation (which is not a pure mode eigenfunction) 
quickly reaches the typical form of the eigenfunctions, and be- 
gins to grow rapidly on the expected timescale. 

We note the similarity of the field structure at f = 0.01 - 0.02 
with Figs. 3 and 4 in RG, which delineate the structure of the 
most rapidly growing mode. During the linear stage (a few Hall 
times), the unstable modes grow exponentially retaining their 
structure, as expected, until they reach amplitudes large enough 
for non-linear effects to be important. At that point (> 0.02), 
the evolution becomes non-linear, there is a drift towards the in- 
ner boundary, and the global shape of the field is completely 
reorganized. After a stationary configuration is adopted, we en- 
ter the diffusion regime and observe long-lived, damped oscilla- 
tory modes with smaller scale modes that dissipate faster During 
the late evolution (bottornpanel of Fig. [TJ, only the longest 
wavelength modes survive [J The evolution of the field in terms 
of the power spectra is shown in Fig. |2l The fast growth of 
the initial perturbation (fcj = 0.5) and other unstable modes 
(fcf = 1.0, 1.5, ...) at early times is evident. After the unstable- 
mode energy has grown several orders of magnitude, non-linear 
coupling fills the full spectrum by f = 0.03, but the dominant 
mode (besides k^ - 0) is always the k^ - 0.5 mode, which 
has drained more than 10% of the background energy. Collapsed 
onto the k- axis, the spectra only shows the ordinary cascade fill- 
ing short wavelength modes, but this only indicates that the field 
has drifted vertically, been compressed into a narrower region, 
and should not be interpreted as extended turbulence over the 
whole domain. 

In contrast to WH, who never detected the growth of the ini- 
tial perturbation and thus concluded that it is always subsumed 
into the standard Hall cascade, we observed the same unstable 
behavior shown in Figs. [T]and|2]in a number of other models 
with different amplitude and k^ of the initial perturbation. The 
same result is obtained with different grid resolutions. We also 
explored other initial choices for f(z), (e.g. f(z) oc sin(fe.z)) and 
replaced the outer boundary condition with the simpler condition 
Bx - By - Q obtaining qualitatively similar results. What then 
is the origin of this discrepancy? We note that the finite size (L) 
of the grid fixes the lower limit of ky. Periodic boundary condi- 
tions result in numerical truncation of all modes with ky < 2n/L. 
In the models in which we restricted the horizontal grid to the 
interval [0, 27r] (the typical size in WH models), we observed 
that unstable modes grew initially but the poloidal perturbations 
reached only a few percent of the total energy and the toroidal 
component saturated even earlier. Although the linear analysis 
predicted a few unstable modes with ky > 1, in this restricted 
domain that does not allow modes with ky < 1 the non-linear 
coupling together with the truncation of long wavelength modes 
somehow resulted in the artificial damping of unstable modes 
much earlier than expected. However, a simple increase in the 
grid size L by a factor of 2, allowing for extra modes to exist 
{ky = 0.5 but also k^ - 1 .5, 2.5...), immediately changed the pic- 
ture, allowing the perturbations to grow to ~ Bq. WH included 
a large number of modes but with > 1, thus the truncation 
of longer wavelength modes seems to be the reason for the ar- 
tificial damping of the HI. We note that, as the magnetic field 
geometry deviates from its initial shape, either due to the finite 
resistivity (the background field is not affected by the Hall drift) 
or because of the growth of the perturbations, the linear analy- 



' Animations of the long-term simulations are available at 
http : / /www . ua . es/per sonal/ j ose . pons/movi es . html , where 
the evolution is visualized in much more detail than a figure can show. 
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Fig. 2. Evolution of the power spectrum. Spectra are shown at 
times t - 0,0.01,0.02 and 0.03 with asterisks, solid lines, dots, 
and dashes, respectively. The total energy is collapsed onto the 
kx axis (top) and k, axis (bottom). The upper dotted line marks 
the initial power in the k^- = background field, from which the 
energy of the unstable modes is drained. 



of a longitudinal size 70 km and a thickness of 1 km (but 
with reflecting boundary conditions instead of periodic bound- 
ary conditions at the poles as used in this work). This resembles 
one of our models that has the largest grid size and therefore 
for which the HI occurs. By analogy with this toy model, one 
may be tempted to claim that the HI is truly relevant and will 
affect the evolution of the magnetic field in NSs. When we com- 
pare the structures and evolutions of the models shown in this 
paper with those of Pons & Geppert (2007), both the qualitative 
shape of the field and the typical evolution are very similar, de- 
spite important differences in the physical model. Perhaps the 
octupolar structures that survive to la te times in the realistic NS 
crust simulations of [Pons Geppert) (2007) originate from the 
HI. However, caution is required in accepting this interpretation. 
Given the sensitivity of the results to details such as the size of 
the grid (which allows only a finite number of unstable modes) 
and the existence of non-linear couplings between modes, it be- 
comes necessary to perform multidimensional simulations in a 
spherical shell with a realistic crust model before we can make 
any claim about the relevance of the HI in NSs. In addition, toy 
models that have constant density and cubic boxes with periodic 
boundary conditions, such as those studied in WH, cannot be 
used to argue against the relevance of the HI in NSs, which is a 
global, more complex problem, that also depends on the true di- 
mensions and the nature of the physical boundaries. As we have 
proven, the numerical truncation of some modes due to the com- 
putational domains being too small can suppress the instability. 
An extension of our present code to spherical coordinates and to 
contain a realistic NS model is under development. Very prelim- 
inary results seem to indicate that an active Hall stage exists, but 
only a more detailed analysis will provide a definitive answer 
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sis is not longer valid. As the field changes, the /:^-window of 
instability may shift to lower values, and in that case the insta- 
bility will be suppressed unless modes with longer wavelengths 
are permitted. Another important feature confirming the validity 
of our results is that the linear stage always proceeds similarly, 
irrespectively of the size of the initial perturbation (if models are 
compared with the appropriate scaling factor). 

4. Conclusions 

We have proven that there are cases in which the HI dominates 
the early evolution of the magnetic field and completely reorga- 
nizes the field structure in a few Hall times. It is not always sub- 
sumed by the Hall cascade, as argued by other authors. We have 
found that the finite size of the computational domain, and the 
use of different boundary conditions significantly aff'ect the out- 
come. Since the HI is not entirely a scale-free problem and, when 
a slab geometry is considered, the relative sizes of the length and 
thickness of the slab make a difference, our findings raise new 
and interesting questions. In simple terms, a thin slab is easier 
subject to the HI than a thicker one but, what happens in a real 
NS crust ? 

The crust of a NS is not an infinite slab, but a finite-size, 
curved domain. To gain insight into realistic scenarios, one can 
mimic the axially symmetric case with a plane parallel model 
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Fig. 1. Evolution of the initial configuration defined by Eq. |2]with a - 0.5 perturbation at four different times. The left panels 
show transverse cuts of the By (dashes) and B, (solid line) components at x = 0, while the right panels show the field structure. 
Arrows show the B^, B^-components and the color scale the By component (red/black positive, yellow/white negative). Note the 
different scales in the left panels. 



